St. Jude OS data CNV analysis

St. Jude OS data CNV analysis

Goals:

We want to characterize the CNV data from St. Jude OS samples in order to determine if the CNVs are consistent across samples within a patient. This will tell us if these tumors are genomically stable or if there are ongoing rearrangements. We also want to determine if there are sub-clones within a single patient which would be revealed by distinct CNV patterns within a patient.

Samples:

For this analysis we got access to many St. Jude OS samples. We also included our single-cell sequencing data which we treated as bulk for the purposes of these analyses.

Workflow

Load libraries

library(tidyverse)

## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'pillar'

## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'hms'

## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──

## ✔ ggplot2 3.3.4     ✔ purrr   0.3.4
## ✔ tibble  3.1.7     ✔ dplyr   1.0.7
## ✔ tidyr   1.1.3     ✔ stringr 1.4.0
## ✔ readr   2.0.1     ✔ forcats 0.5.1

## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

theme_set(theme_bw())
theme_update(plot.title = element_text(hjust = 0.5))
library(gt)

Make up directory structure

for directoryName in \
  misc \
  slurmOut \
  input\split \
  output \
  output/aligned \
  output/figures \
  output/counts \
  output/varscan/copynumber/ \
  output/vcfs \
  output/vcfs/realignSummary \
  output/vcfs/p53Char \
  output/vcfs/p53Char/clinVar
do
  if [ ! -d ${directoryName} ]
  then
    mkdir -p ${directoryName}
  fi 
done

Prepare files that contain links to DNAnexus data

Had to do this by hand due to the links expiring after 24 hours

We had a lot of trouble downloading large bam files from DNAnexus. The downloads frequently got interupted. To get around this, I split the files in DNAnexus into smaller chunks to make interuptions less likely. I also wrote the code below to automatically resume the downloads if the download was interupted. I had the script check the md5sums to be sure that the file wasn’t corrupted during download. The original file is then re-generated using cat and md5sum checked again.

Don’t forget to dos2unix the .txt files before running!

perl ../scripts/prepDownloadFile.pl \
  --linkFile misc/DNAnexus_export_urls-20210926-103118.txt \
  --md5File misc/allMd5.txt \
    > misc/batch1DownloadFile.txt

perl ../scripts/prepDownloadFile.pl \
  --linkFile misc/DNAnexus_export_urls-20210927-085055.txt \
  --md5File misc/allMd5.txt \
    > misc/batch2DownloadFile.txt

perl ../scripts/prepDownloadFile.pl \
  --linkFile misc/DNAnexus_export_urls-20210928-082401.txt \
  --md5File misc/allMd5.txt \
    > misc/batch3DownloadFile.txt

perl ../scripts/prepDownloadFile.pl \
  --linkFile misc/DNAnexus_export_urls-20210929-083226.txt \
  --md5File misc/allMd5.txt \
    > misc/batch4DownloadFile.txt

perl ../scripts/prepDownloadFile.pl \
  --linkFile misc/DNAnexus_export_urls-20210930-084226.txt \
  --md5File misc/allMd5.txt \
    > misc/batch5DownloadFile.txt
    
perl ../scripts/prepDownloadFile.pl \
  --linkFile misc/DNAnexus_export_urls-20211001-084849.txt \
  --md5File misc/allMd5.txt \
    > misc/batch6DownloadFile.txt

Download split files from DNA nexus

I had to run these one by one manually due to how the data are organized

sbatch sbatchCmds/forceDownload_1.sh
sbatch sbatchCmds/forceDownload_2.sh
sbatch sbatchCmds/forceDownload_3.sh
sbatch sbatchCmds/forceDownload_4.sh
sbatch sbatchCmds/forceDownload_5.sh
sbatch sbatchCmds/forceDownload_6.sh

Downloaded the md5 data for the raw bam files

cat misc/rawMd5s/rawBamMd5_md5sum_job-G5FQ* > misc/rawBamMd5s.txt

Recombine the split files and check md5sums

sbatch sbatchCmds/catSplitFiles.sh

Make up symbolic links to previously analyzed data so we can analyze those as well

for original in S0113 S0114 S0115 S0116 S0126 S0127 S0128
do
    ln -s \
        /home/gdrobertslab/lab/Counts/${original}/possorted_bam.bam \
        input/${original}.bam
done

# Some of the samples were already in the input folder and don't need symLinks
for original in \
        SJOS046149_R1 SJOS046149_R2 SJOS046149_X1 SJOS031478_D1 \
        SJOS031478_D2 SJOS031478_D3
do
    ln -s \
        ~/analyses/roberts/stjude/rawData/${original}.WholeGenome.bam \
        input/${original}.WholeGenome.bam
done

Make up a list of all chromosomes in the hg38 genome

Kicked out the unknown and random chromosomes

cat /reference/homo_sapiens/hg38/ucsc_assembly/illumina_download/Sequence/WholeGenomeFasta/genome.fa.fai | \
  grep -v "chrUn\|random\|EBV" \
    > misc/chrList.txt

Realign the bam files using a mixed human/mouse reference

if [ ! -d /gpfs0/scratch/mvc002/roberts/stjude/realigned ]
then
    mkdir -p /gpfs0/scratch/mvc002/roberts/stjude/realigned
fi

sbatch sbatchCmds/realign.sh

Call CNVs using varscan

sbatch sbatchCmds/mpileupVarscan.sh

Refine CNV calls using varscan copycaller

sbatch sbatchCmds/varScanCopycaller.sh

Plot the distribution of CNVs found by copynumber

pdf(file = "output/figures/CNVdistribution_copynumber.pdf")
for (sample_name in list.dirs(path = "output/varscan/copynumber",
                              full.names = FALSE,
                              recursive = FALSE)) {
    print(read_tsv(list.files(path = paste("output/varscan/copynumber/",
                                           sample_name,
                                           sep = ""),
                             pattern = "*.copynumber",
                             full.names = TRUE),
                   show_col_types = FALSE) %>%
        sample_n(1000000) %>%
        select(log2_ratio) %>%
        mutate(ratio = 2^log2_ratio) %>%
        ggplot(aes(x = ratio)) +
        geom_histogram(bins = 300) +
        geom_vline(xintercept = 1) +
        scale_x_log10() +
        geom_vline(xintercept = 2) +
        geom_vline(xintercept = 0.5) +
        ggtitle(sample_name))
}
dev.off()

## png 
##   2

Calculate the distance between samples based on CNVs

This approach has the issue that normalization isn’t perfect across the samples due to the extensive number of genomic changes. This artificially increases the distance between similar samples. This likely will not be used in the paper.

perl scripts/varscanToMatrix.pl \
    -v \
    --fileList "output/varscan/copycaller/S*/*chr*.txt" \
        > output/varscan/varscan_distMatrix.txt

perl scripts/varscanToMatrix_threshhold.pl \
    -v \
    --fileList "output/varscan/copycaller/S011*/*chr*.txt" \
        > output/varscan/varscan_distMatrix_threshhold.txt

Combine the varscan outputs into one file

sbatch sbatchCmds/combineVarscan.sh

Make plots of the CNV calls across all the samples

large_bin_size <- 100000

combined_data <- read_tsv("output/varscan/combined_calls.txt",
                          show_col_types = FALSE) %>%
    mutate(chr_num = str_remove(chr, "chr") %>%
                as.numeric()) %>%
    arrange(chr_num, bin) %>%
    mutate(order = seq_len(nrow(.)),
           larger_bin = floor(bin / large_bin_size) * large_bin_size) %>%
    pivot_longer(cols = c(-chr, -bin, -order, -chr_num, -larger_bin),
                 names_to = "sample",
                 values_to = "log_ratio")

patient_key <- data.frame(
    patient = unique(combined_data$sample) %>%
        str_remove("_.+") %>%
        str_replace("S0126", "SJOS030645") %>%
        str_replace("S0127", "SJOS030645") %>%
        str_replace("S0128", "SJOS030645") %>%
        str_replace("S0114", "SJOS031478") %>%
        str_replace("S0116", "SJOS031478") %>%
        str_replace("S0113", "SJOS046149") %>%
        str_replace("S0115", "SJOS046149"),
    sample_type = unique(combined_data$sample) %>%
        str_replace("^S0.+", "Xenograft") %>%
        str_replace(".+_X.+", "Xenograft") %>%
        str_replace(".+_D.+", "Diagnosis") %>%
        str_replace(".+_M.+", "Metastasis") %>%
        str_replace(".+_R.+", "Relapse"),
    sample = unique(combined_data$sample))

combined_data <- full_join(combined_data, patient_key)

## Joining, by = "sample"

chr_cols <- sample(rainbow(length(unique(combined_data$chr))),
                   length(unique(combined_data$chr)),
                   replace = FALSE)

for (patient_name in unique(patient_key$patient)) {
    plot_name <- combined_data %>%
        filter(patient == patient_name) %>%
        group_by(chr, sample, larger_bin) %>%
        summarize(log_ratio = median(log_ratio),
                order = min(order),
                .groups = "drop") %>%
        ggplot(aes(x = order,
                    y = log_ratio,
                    color = chr)) +
        geom_hline(yintercept = 0, color = "black") +
        geom_point(alpha = 0.5,
                    size = 0.5) +
        scale_color_manual(values = chr_cols) +
        labs(y = "Copy number log ratio",
             x = "") +
        facet_wrap(~ sample, ncol = 1)

    ggsave(paste0("output/figures/varscan",
                 patient_name,
                 ".png"),
           plot = plot_name,
           width = 20,
           height = 8)
}

## Take a look at regions commonly amplified/deleted across samples
unique_clones <- c("SJOS001101_M4",
                   "SJOS001111_M1",
                   "SJOS001116_X2",
                   "SJOS001121_X2",
                   "SJOS001126_X1",
                   "SJOS013768_X1",
                   "SJOS016016_X1",
                   "SJOS030101_R1",
                   "SJOS030645_D2",
                   "SJOS031478_D3",
                   "SJOS046149_R2",
                   "SJOS001105_R1")

combined_data %>%
    filter(sample %in% unique_clones) %>%
    group_by(chr, bin) %>%
    summarize(mean_ratio = median(log_ratio),
              order = min(order),
              .groups = "drop") %>%
    ggplot(aes(x = bin,
               y = mean_ratio,
               color = chr)) +
    geom_hline(yintercept = 0, color = "black") +
    geom_point(alpha = 0.5,
               size = 0.5) +
    labs(x = "",
         y = "Mean copy number ratio") +
    facet_wrap(~ chr, scales = "free_x") +
    theme(legend.position = "none")

ggsave("output/figures/clonal_mean_CNV.png",
       width = 15,
       height = 10)

high_chr3 <- combined_data %>%
    filter(sample %in% unique_clones) %>%
    group_by(chr, bin) %>%
    summarize(mean_ratio = median(log_ratio),
              order = min(order),
              .groups = "drop") %>%
    filter(chr == "chr3" & mean_ratio < -0.5)

Look at correlation of CNVs between samples

cnv_cor <- read_tsv("output/varscan/combined_calls.txt",
                     col_names = TRUE,
                     show_col_types = FALSE) %>%
    mutate(chr_bin = paste(chr, bin, sep = "_"), .keep = "unused") %>%
    column_to_rownames("chr_bin") %>%
    cor()

patient_key <- data.frame(patient = rownames(cnv_cor) %>%
                                str_remove("_.+") %>%
                                str_replace("S0126", "SJOS030645") %>%
                                str_replace("S0127", "SJOS030645") %>%
                                str_replace("S0128", "SJOS030645") %>%
                                str_replace("S0114", "SJOS031478") %>%
                                str_replace("S0116", "SJOS031478") %>%
                                str_replace("S0113", "SJOS046149") %>%
                                str_replace("S0115", "SJOS046149"),
                          sample_type = rownames(cnv_cor) %>%
                                str_replace("^S0.+", "Xenograft") %>%
                                str_replace(".+_X.+", "Xenograft") %>%
                                str_replace(".+_D.+", "Diagnosis") %>%
                                str_replace(".+_M.+", "Metastasis") %>%
                                str_replace(".+_R.+", "Relapse"),
                          sample_2 = rownames(cnv_cor)) %>%
    column_to_rownames("sample_2")

long_cor <-
    cnv_cor %>%
    as.data.frame() %>%
    rownames_to_column("sample_1") %>%
    pivot_longer(-sample_1,
                 names_to = "sample_2",
                 values_to = "cor") %>%
    left_join(patient_key %>%
                rownames_to_column("sample_1") %>%
                rename(patient_1 = patient) %>%
                select(-sample_type)) %>%
    left_join(patient_key %>%
                rownames_to_column("sample_2") %>%
                rename(patient_2 = patient) %>%
                select(-sample_type)) %>%
    mutate(within = patient_1 == patient_2,
           sample_1 = fct_reorder(sample_1, patient_1),
           sample_2 = fct_reorder(sample_2, patient_2)) %>%
    filter(sample_1 != sample_2) %>%
    arrange(patient_1)

## Joining, by = "sample_1"Joining, by = "sample_2"

## Try to id clones
clone_table <- tibble(sample_1 = character(),
                      clone_num = factor())
claimed_list <- character()
cor_cutoff <- 0.4
clone_num <- 1

for (sample_name in rownames(cnv_cor)) {
    if (!sample_name %in% claimed_list) {
        new_claimed <-
            long_cor %>%
            filter(sample_1 == sample_name &
                    cor >= cor_cutoff &
                    within == TRUE) %>%
            select(sample_2, cor) %>%
            mutate(clone_num = as.factor(clone_num)) %>%
            bind_rows(tibble(clone_num = as.factor(clone_num),
                             sample_2 = sample_name,
                             cor = 1))
        clone_table <-
            rbind(clone_table, new_claimed)
        claimed_list <- c(claimed_list, sample_name, new_claimed$sample_2)
        clone_num <- clone_num + 1
    }
}


clone_table <-
    clone_table %>%
    group_by(sample_2) %>%
    arrange(cor * -1, .by_group = TRUE) %>%
    slice_head(n = 1) %>%
    ungroup() %>%
    select(-cor)

long_cor <- full_join(long_cor, clone_table, by = "sample_2")

# Mean within patient correlation
long_cor %>%
    filter(within == TRUE) %>%
    pull(cor) %>%
    summary()

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.007455 0.362223 0.509682 0.503482 0.717656 0.921232

# # Mean correlation within a clone
# long_cor %>%
#     select(sample_2, patient_2, cor, clone_num) %>%
#     group_by(patient_2, clone_num) %>%
#     summarize(cor = mean(cor)) %>%
#     pull(cor) %>%
#     summary()

# Get number of patients with >1 clone
long_cor %>%
    select(patient_2, clone_num) %>%
    arrange(patient_2) %>%
    distinct() %>%
    pull(patient_2) %>%
    table()

## .
## SJOS001101 SJOS001105 SJOS001111 SJOS001116 SJOS001121 SJOS001126 SJOS013768 
##          1          1          1          1          1          1          2 
## SJOS016016 SJOS030101 SJOS030645 SJOS031478 SJOS046149 
##          1          2          2          2          2

# Mean within clone correlation
long_cor %>%
    full_join(clone_table %>%
                rename(sample_1 = sample_2,
                       clone_2 = clone_num)) %>%
    filter(clone_num == clone_2) %>%
    pull(cor) %>%
    summary()

## Joining, by = "sample_1"

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4007  0.5312  0.6807  0.6707  0.7850  0.9212

# Mean correlatin within a patient between divergent clones
long_cor %>%
    full_join(clone_table %>%
                rename(sample_1 = sample_2,
                       clone_2 = clone_num)) %>%
    filter(patient_1 == patient_2 &
           clone_num != clone_2) %>%
    pull(cor) %>%
    summary()

## Joining, by = "sample_1"

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.007455 0.123653 0.307508 0.282802 0.409872 0.584116

# Mean correlation between patients
long_cor %>%
    full_join(clone_table %>%
                rename(sample_1 = sample_2,
                       clone_2 = clone_num)) %>%
    filter(patient_1 != patient_2) %>%
    pull(cor) %>%
    summary()

## Joining, by = "sample_1"

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.1363  0.1165  0.1823  0.1781  0.2438  0.4527

# Correlation between xenografts and original samples within a single clone and patient
long_cor %>%
    full_join(clone_table %>%
                rename(sample_1 = sample_2,
                       clone_2 = clone_num)) %>%
    left_join(patient_key %>%
                rownames_to_column("sample_1") %>%
                rename(type_1 = sample_type) %>%
                select(-patient)) %>%
    left_join(patient_key %>%
                rownames_to_column("sample_2") %>%
                rename(type_2 = sample_type) %>%
                select(-patient)) %>%
    filter(patient_1 == patient_2 &
           clone_num == clone_2 &
           type_1 == "Xenograft" &
           type_2 != "Xenograft") %>%
    select(-within, -patient_1, -patient_2) %>%
    pull(cor) %>%
    summary()

## Joining, by = "sample_1"Joining, by = "sample_1"Joining, by = "sample_2"

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4710  0.5622  0.7224  0.7019  0.8295  0.9212

ggplot(long_cor, aes(x = cor)) +
    geom_histogram(bins = 200) +
    facet_wrap(~ within, ncol = 1)

# correlation of about 0.4 seems to be the upper bound of non-related samples

n_clones <- max(as.numeric(long_cor$clone_num))

text_col <-
    long_cor %>%
    full_join(tibble(clone_num = seq_len(n_clones) %>%
                        as.factor(),
                    text_col = rainbow(n = n_clones) %>%
                        sample(size = n_clones))) %>%
    mutate(group = paste(patient_2, sample_2)) %>%
    select(group, text_col) %>%
    distinct() %>%
    arrange(group) %>%
    pull(text_col)

## Joining, by = "clone_num"

ggplot(long_cor, aes(x = paste(patient_1, sample_1),
                     y = paste(patient_2, sample_2),
                     fill = cor)) +
    geom_tile() +
    scale_fill_gradient2(low = "blue",
                         mid = "white",
                         high = "red",
                         midpoint = 0.5) +
    theme(axis.text.y = element_text(colour = text_col))

## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

patient_key <- patient_key %>%
    rownames_to_column("sample_2") %>%
    full_join(clone_table) %>%
    column_to_rownames("sample_2") %>%
    rename(clone = clone_num,
           `sample type` = sample_type) %>%
    relocate(`sample type`, .after = last_col())# %>%

## Joining, by = "sample_2"

#    select(-`sample type`)

col_patient <- list(patient = c("#a6cee3",
                                "#1f78b4",
                                "#b2df8a",
                                "#33a02c",
                                "#fb9a99",
                                "#e31a1c",
                                "#fdbf6f",
                                "#ff7f00",
                                "#cab2d6",
                                "#6a3d9a",
                                "#ffff99",
                                "#b15928"),
                    clone = c("#96FF00",
                              "#00FFD2",
                              "#FF00B4",
                              "#2b754e",
                              "#3CFF00",
                              "#FF5A00",
                              "#F000FF",
                              "#9600FF",
                              "#3C00FF",
                              "#F0FF00",
                              "#FFB400",
                              "#0078FF",
                              "#001EFF",
                              "#00FF1E",
                              "#FF0000",
                              "#00D2FF",
                              "#FF005A"),
                    `sample type` = c("#a6cee3",
                                    "#1f78b4",
                                    "#b2df8a",
                                    "#33a02c"))

names(col_patient$patient) <- patient_key$patient %>%
    as.factor() %>%
    levels()

names(col_patient$`sample type`) <- patient_key$`sample type` %>%
    as_factor() %>%
    levels()

names(col_patient$clone) <- levels(clone_table$clone_num)

png("output/figures/cnv_correlation_heatmap.png",
    width = 3500,
    height = 3000,
    res = 300)
pheatmap::pheatmap(cnv_cor,
                   annotation_row = patient_key,
                   annotation_col = patient_key,
                   annotation_colors = col_patient)
dev.off()

## png 
##   3

SNP calling

Call snps using bcftools mpileup and filter SNPs

Filter out indels and SNPs with qual < 20 and depth < 20

# Make up a list of chromosomes to analyze for use by parallel
grep "^>" \
    /reference/homo_sapiens/hg38/ucsc_assembly/illumina_download/Sequence/WholeGenomeFasta/genome.fa \
    | grep -v "random\|chrUn\|chrEBV" \
    | perl -pe 's/>//' \
        > misc/chrList.txt

sbatch sbatchCmds/callVariants.sh

sbatch sbatchCmds/indexVcfs.sh

Calculate distance between samples based on SNP calls

module load GCC/9.3.0 \
            GCCcore/9.3.0 \
            BCFtools/1.11 \
            SAMtools/1.15

perl scripts/vcfToMatrix_parallel.pl \
    --fileList "output/vcfs/S*vcf.gz" \
    --threads 20 \
        > output/vcfs/snpDistanceMatrix.txt

Check P53 for mutations using SnpEff

sbatch sbatchCmds/TP53MutChar.sh

wget https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar_20220416.vcf.gz
mv clinvar_20220416.vcf.gz misc/
zcat misc/clinvar_20220416.vcf.gz \
    | awk '$1 ~ /^#/ || \
           $1 == 17 && \
           $2 >= 7600000 && \
           $2 <= 7690000 \
           {print}' \
    | perl -pe 's/^17/chr17/' \
    | gzip \
    > misc/clinvar_TP53.vcf.gz

Merge clinvar data with TP53Muts.txt

combineClinvar.sh

filter_list <- function(x) {
    str_split(x, ",") %>%
    unlist() %>%
    str_subset("TP53") %>%
    str_c(collapse = ",")
}

file_list <- list.files("output/vcfs/p53Char/clinVar",
                        pattern = ".+_G.+",
                        full.names = TRUE)

combined_data <- tibble(sample = character())

for (file_name in file_list) {
    combined_data <-
        read_tsv(file_name,
                 col_names = TRUE,
                 col_types = cols(.default = "c")) %>%
        filter(grepl("TP53", ANN)) %>%
        rowwise() %>%
        mutate(ANN = filter_list(ANN),
               sample = str_match(file_name, ".+/(.+).txt")[2]) %>%
        filter(ANN != "") %>%
        bind_rows(combined_data) %>%
        relocate(sample, CLNSIG, ANN)
}

write_tsv(combined_data,
          file = "output/vcfs/p53Char/p53Char_combined_germline.txt")

perl scripts/cmpVarscanCalls.pl \
    --binSize 100 \
    --inputFiles "output/varscan/copycaller/*/*chr17.txt" \
        > output/varscan/combined_calls_chr17.txt

## 
## Printing out results

Plot P53 CNV data

P53_loc <- list(min = 7661779,
                max = 7687538)

combined_data <- read_tsv("output/varscan/combined_calls_chr17.txt",
                          show_col_types = FALSE) %>%
    mutate(chr_num = str_remove(chr, "chr") %>%
                as.numeric()) %>%
    filter(chr_num == 17 & bin >= 7550000 & bin <= 7800000) %>%
    arrange(chr_num, bin) %>%
    pivot_longer(cols = c(-chr, -bin, -chr_num),
                 names_to = "sample",
                 values_to = "log_ratio")

patient_key <- data.frame(
    patient = unique(combined_data$sample) %>%
        str_remove("_.+") %>%
        str_replace("S0126", "SJOS030645") %>%
        str_replace("S0127", "SJOS030645") %>%
        str_replace("S0128", "SJOS030645") %>%
        str_replace("S0114", "SJOS031478") %>%
        str_replace("S0116", "SJOS031478") %>%
        str_replace("S0113", "SJOS046149") %>%
        str_replace("S0115", "SJOS046149"),
    sample_type = unique(combined_data$sample) %>%
        str_replace("^S0.+", "Xenograft") %>%
        str_replace(".+_X.+", "Xenograft") %>%
        str_replace(".+_D.+", "Diagnosis") %>%
        str_replace(".+_M.+", "Metastasis") %>%
        str_replace(".+_R.+", "Relapse"),
    sample = unique(combined_data$sample))

combined_data <- full_join(combined_data, patient_key) %>%
    mutate(amplification = if_else(log_ratio > 0.5,
                                   "amplification",
                                   if_else(log_ratio < -0.5,
                                           "deletion",
                                           "normal")))

## Joining, by = "sample"

ggplot(combined_data,
       aes(x = bin,
           y  = sample,
           fill = log_ratio)) +
    geom_tile() +
    scale_fill_gradient2(low = "#8e0152",
                         mid = "white",
                         high = "#276419") +
    facet_grid(patient ~ .,
               scales = "free_y",
               space = "free") +
    theme(strip.text.y.right = element_text(angle = 0)) +
    geom_vline(xintercept = P53_loc$min,
               linetype = "dashed",
               colour = "red") +
    geom_vline(xintercept = P53_loc$max,
               linetype = "dashed",
               colour = "red") +
    labs(x = "Chromosome 17 position",
         y = "",
         title = "TP53 copy numbers")

ggsave(file = "output/figures/p53Cnv_chr17.png",
       width = 10,
       height = 10)

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux 8.1 (Ootpa)
## 
## Matrix products: default
## BLAS:   /gpfs0/export/apps/opt/R/4.1.0-foss-2020a/lib64/R/lib/libRblas.so
## LAPACK: /gpfs0/export/apps/opt/R/4.1.0-foss-2020a/lib64/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gt_0.3.0        forcats_0.5.1   stringr_1.4.0   dplyr_1.0.7    
##  [5] purrr_0.3.4     readr_2.0.1     tidyr_1.1.3     tibble_3.1.7   
##  [9] ggplot2_3.3.4   tidyverse_1.3.1 rmarkdown_2.9  
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.1  xfun_0.24         bslib_0.2.5.1     haven_2.4.1      
##  [5] colorspace_2.0-1  vctrs_0.3.8       generics_0.1.0    htmltools_0.5.1.1
##  [9] yaml_2.2.1        utf8_1.2.1        rlang_1.0.2       jquerylib_0.1.4  
## [13] pillar_1.7.0      glue_1.4.2        withr_2.4.2       DBI_1.1.1        
## [17] dbplyr_2.1.1      modelr_0.1.8      readxl_1.3.1      lifecycle_1.0.0  
## [21] munsell_0.5.0     gtable_0.3.0      cellranger_1.1.0  rvest_1.0.0      
## [25] evaluate_0.14     knitr_1.33        tzdb_0.1.2        ps_1.6.0         
## [29] fansi_0.5.0       broom_0.7.7       Rcpp_1.0.7        backports_1.2.1  
## [33] scales_1.1.1      jsonlite_1.7.2    fs_1.5.0          hms_1.1.0        
## [37] digest_0.6.27     stringi_1.6.2     grid_4.1.0        cli_2.5.0        
## [41] tools_4.1.0       magrittr_2.0.1    sass_0.4.0        crayon_1.4.1     
## [45] pkgconfig_2.0.3   ellipsis_0.3.2    xml2_1.3.2        reprex_2.0.0     
## [49] lubridate_1.7.10  rstudioapi_0.13   assertthat_0.2.1  httr_1.4.2       
## [53] R6_2.5.0          compiler_4.1.0